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RECONSTRUCTION OF SIGNALS FROM MAGNITUDES OF 
REDUNDANT REPRESENTATIONS: THE COMPLEX CASE 

RADU BALAN 



Abstract. This paper is concerned with the question of reconstructing a vector in a finite- 
dimensional complex Hilbert space when only the magnitudes of the coefficients of the vector 
under a redundant linear map are known. We present new invertibility results as well an 
iterative algorithm that finds the least-square solution and is robust in the presence of noise. 
.^ We analyze its numerical performance by comparing it to the Cramer-Rao lower bound. 

^D 

<T^ 1. Introduction 

This paper is concerned with the question of reconstructing a vector a; in a finite-dimensional 
complex Hilbert space H when only the magnitudes of the coefficients of the vector under 
a redundant linear map are known. 

Specifically our problem is to reconstruct x G H up to a global phase factor from the 

magnitudes {\{x, fk)\ , 1 < k < m} where {/i, • • • , f m } is a frame (complete system) for H. 

This paper is the companion paper to the arXiv submission [9]. 

Q\ A previous paper [3] described the importance of this problem to signal processing, in 

particular to the analysis of speech. Of particular interest is the case when the coefficients 

are obtained from a Windowed Fourier Transform (also known as Short-Time Fourier Trans- 

■^ form), or an Undecimated Wavelet Transform (in audio and image signal processing). While 

[3] presents some necessary and sufficient conditions for reconstruction, the general problem 

of finding fast/efficient algorithms is still open. In [8] we describe one solution in the case 

of STFT coefficients. The authors of [57] propose a novel algorithm adapted to compactly 

supported signals and FFT transforms that uses individual signal spectral powers and of two 

rS additional interferences between signals. 

For vectors in real Hilbert spaces, the reconstruction problem is easily shown to be equiv- 
alent to a combinatorial problem. In [5] this problem is further proved to be equivalent to a 
(nonconvex) optimization problem. 

A different approach (which we called the "algebraic approach") was proposed in [7J. While 
it applies to both real and complex cases, noisless and noisy cases, the approach requires 
solving a linear system of size exponentially in space dimension. The algebraic approach 
mentioned earlier generalizes the approach in [6] where reconstruction is performed with 
complexity 0(n 2 ) (plus computation of the principal eigenvector for a matrix of size n) but 
requiring m = 0(n 2 ) frame vectors. 

Recently the authors of [SU] developed a convex optimization algorithm (PhaseLift) and 
proved its ability to perform exact reconstruction in the absence of noise, as well as its 

l 
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stablity under noise conditions. In a separate paper, [21], the authors developed further 
a similar algorithm in the case of windowed DFT transforms. The original requirement of 
m = 0(nlog(n)) vectors has been relaxed to m = 0(n) in [22]- Similar convex optimization 
solutions have been proposed by the authors of [29] and [61]. Additionally, [36J studied 
the duality gap in this approach and obtained a necessary and sufficient condition for the 
existence of a dual certificate. 

While writing this paper, we have become aware of [UJ where the authors obtained similar 
results to the injectivity criteria presented here, as well as to the Cramer-Rao Lower Bound 
derived in this paper. We will comment more later in the paper. We also acknowledge the 
paper [30] where some stability bounds have been obtained in the real case. The real case 
of the stabilty bounds obtained here are presented in a separate paper [TO] together with 
additional results for the real case. 

The organization of the paper is as follows. In section [2] we define the problem explicitely. 
In section [3] we describe new analysis results. Specifically we analysis in more details spaces 
of symmetric operators of constrained signature, and then we shows how they are related to 
the phaseless recovery problem. Our results are canonical, meaning that are independent to a 
particular choice of basis. In section [4] we establish two robustness results: bi-Lipschitzianity 
of the nonlinear analysis map, and the Cramer- Rao lower bound (CRLB). In section p^ we 
present a new reconstruction algorithm based on the Least-Square method. We also obtain 
robustness bounds to noise. Its performance is analyzed in section [6] and is compared to the 
CRLB. Section Section [7] contains conclusions and is followed by references. 

2. Background 

Let H be an n-dimensional complex Hilbert space (such as C n or C PlXp2 ) with scalar 
product (, ) linear in the first term and antilinear in the second term and a conjugation 
c : H —7- H (see e.g. [35] : conjugation is an antilinear transformation thatt squares to 
identity). Let F = {/i, • • • , f m } be a spanning set of m vectors in H. In finite dimension (as 
it is the case here) such a set forms a frame. In the infinite dimensional case, the concept 
of frame involves a stronger property than completness (see for instance [23])- We review 
additional terminology and properties which remain still true in the infinite dimensional 
setting. The set F is frame if and only if there are two positive constants < A < B < oo 
(called frame bounds) so that 

m 

A\\x\\ 2 < J2\(xJ k )\ 2 < B\\x\\ 2 
k=i 

When we can choose A = B the frame is said tight. For A = B = 1 the frame is called 
Parseval. A set of vectors F of the n-dimensional Hilbert space H is said to have full spark if 
any subset of n vectors is linearly independent (see [2] for a full discussion of such frames). 

2.1. Problem Definition and Notations. For a vector x G H, the collection of coefficients 
{{x, fk) , 1 < k < m} represents the analysis of vector x given by the frame F. In H we 



RECONSTRUCTION OF SIGNALS FROM MAGNITUDES OF REDUNDANT REPRESENTATIONS: THE COMPLEX Cj 

consider the following equivalence relation: 

(2.1) x,y G H , x ~ y iff y = zx for some scalar z with \z\ = 1 

Note z = e l{p for some real number if. Let H = H/ ~ be the set of classes of equivalence 
induced by this relation. Thus x = {e ia x, < a < 2n}. The analysis map induces the 
following nonlinear map 

(2.2) a:#->(R + ) m , a(x) = (\(x, f k )\ 2 )i< k < m 

where 1R + = {x , x G IR , x > 0} is the set of nonnegative real numbers. In [5] we 
studied when the nonlinear map a is injective, mostly in the real case, and we provided 
some necessary conditions of injectivity in the complex case. We review these results below. 
In this paper we obtain additional injectivity results. We then concentrate on the additive 
white Gaussian noise model 

(2.3) y = a(x) + v , v ~ N(0, a 2 ) 

We describe an algorithm (the Iterative Regularized Least-Square algorithm) to solve the 
estimation problem. We prove some convergence results and we study its performance in 
the noisy case. We shall derive the Cramer-Rao Lower Bound (CRLB) for this model and 
compare the algorithm performance to this bound. 

We describe several objects that will be used in this paper. 

The set B(H) denotes the set of linear operators on H. In B(H) for any 1 < p < oo, the 
p-norm of T e B(H) denoted ||T|| is given by the p-norm of its set of singular eigenvalues. 
In particular for p — 1, the 1-norm ||T'|| 1 is called the nuclear norm of T; for p = 2, the 
2-norm ||T|| 2 = y/tr(T*T) is the Frobenius norm of T; for p = oo, oo-norm HT^ is the 
same as the operator norm of T on H simply denoted ||T||. 

For each frame vector fk we denote by F& its associated rank-1 operator 

(2.4) F k :H^H, F k (x) = (x, f k )f k 

In general the associated rank-1 operator to a vector x G H is the operator X : H — >■ H, 
X(v) = (v,x)x. Note this is a slight misnomer: X has at most rank one. Specifically, X has 
rank one if and only if x ^ 0; otherwise X has rank zero. One can think of X as simply the 
n x n matrix xx*, where x* denotes the transpose and conjugate of x. However, as much as 
possible, we plan to use canonical notations (that is, basis independent). 

For any two vectors u, v G H we define their symmetric outer product denoted [w, v\ by 

(2.5) lu,vj:H^H, [«,v](x) = -((x, u)v + (x, v)u) 

One can think of \u,v\ as the n x n symmetric matrix \{yu* + uv*). Note the rank-1 
associated to a vector x can be written as [a;, a;]. Thus F k = lf k ,f k J. Note also [u,vj is 
IR-bilinear but it is not (C-(bi)linear. Furthermore [ix, vj = lv,uj. 

The nonlinear map a induces a linear map A on the set B{H) of bounded operators on 
H: 

(2.6) A : B{H) -)> C m , {A{T)) k = (Tf k , f k ) = tr{TF k } 
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In general, for two operators T,S G B(H) we denote {T,S) B , H \ = tr{TS*} their Hilbert- 
Schmidt scalar product, where S* is the adjoint of S and tr denotes the trace. Note the scalar 
product is basis independent, but it depends of the underlying Hilbert space structure. Thus 
(A(T))k = (T, F k ) B , H y Also note a(x) = A(X) where X = (•, x)x is the rank-1 associated to 
x. When no danger of confusion we drop the index B(H) from the scalar product notation. 
Let Sym(H) denote set of self-adjoint operators on H, Sym(H) = {T G B(H) , T* = T}. 
We denote by S p ' q or S p ' q (H), the set of self-adjoint operators on H that have at most p 
positive eigenvalues and at most q negative eigenvalues: 

S™ ={T G Sym(H) , Sp(T) = {A 1; • • • , A n } , 

Al > • • • > X q > = \ n -p- q = ■ ■ ■ = A n _g > \ n -q+l > ' ' ' > A n } 

where Sp(T) denotes the spectrum of T (i.e. the set of its eigenvalues). Notice S p,q is not 
a linear space, but instead it is a cone in B(H). This cone property is key in deriving 
robustness and stability bounds later on. 

We denote by A ma;c (T) the largest eigenvalue of T and by A min (T) the smallest eigenvalue 
of T. In particular we are interested in S 1,0 and S 1,1 : 

(2.8) S 1 ' = {Te Sym(H) , rank(T) < 1, A mm (T) = 0} 

(2.9) S 1 ' 1 = {Te Sym(H) , rank(T) < 2, Sp(T) = {X max (T), (n " 2) , A mm (T)}, 

Amax\T) > > \min{T)} 

Note the following obvious inclusions 

(2.10) {0} C S lfi C S 1 ' 1 C Sym(H) , {0} C S ' 1 C 5 1 ' 1 C Sym(H) 

We denote by S p,q the subset of S p ' q of selfadjoint operators that have rank p + q, hence 
exactly p strictly positive eigenvalues and q strictly negative eigenvalues. Thus 

(2.ii) s p > q = s M u s p - l > q U S™" 1 

In particular 

(2.12) S 1 ' 1 = S 1 ' 1 U 5 1 ' U S ' 1 , S l >° = S h0 U {0} 



A more detailed analysis of these sets is presented in subsection |3.1 

Next we describe the realification of the Hilbert space H. To do so canonically we need 
to fix a conjugation c : H — >• H. To the complex Hilbert space H with conjugation c 
we associate its 2n- dimensional real vector space H^ subset of H x H built from vectors 
vr = \{v + c(v )) and vi = ^(v — c(v)) as follows: 

(2-13) H M = | Q(« + c(v)), \.{v- c(v))\ ,veH 

Thus H& is the image of H through the IR-linear map, 

(2.14) j : H -► H x H , j(«) = Q(« + c( v )), i(« - c( v )) 
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Note j is injective and when restricted to j : H — )■ H R it becomes an R- isomorphism. Its 
inverse is given by 

(2.15) r'-Hn^H ,r\u,v)=u + iv 
Let J denote the linear map defined by 

(2.16) J : H x H ^ H x H , J(v, w) = (-w, v) 
Note it is conjugate to the multiplicattion by % in H: 

(2.17) J:H R ^H R , J(S(v))=i(iv) 

Hence H R is J invariant. In H R the induced scalar product is given by 
(2.18) 

(j(u) JH) := {-(v + c(v)), -(w + c{w))) + ( — (v - ic(v)), —(w - c(w))) = real((v, w)) 

We denote by (v,w) M = real((v,w)) the IR-linear inner product on H. Note 

(u, v) = real({u, v)) + i imag({u, v)) = (u, v) R — i{iu, v) R 

= (u,v) R + i(u,iv) M = (j(u),j(v))+i(j(u),Jj(v)) 

If {ei, • • • , e n } is an orthonormal basis in H, then {j(ei), • • • , j(e n ), Jj(ei), • • • , Jj(e n )} is an 
orthonormal basis in H R . Which shoes the real dimension of H R is 2n, dim^ i^jj = In. 
On H x H there are two inner product structures: 

(2.20) ({x,y),(u,v)) c = (x,u) + {y,v) 

(2.21) ((x,y),(u,v)) M = -(((x,y),(u,v)) + ((u,v),(x,y))) = (x,u) m + (y,v) m 

With respect to (•, -) c , H x H is C-vector space of dimension In. With respect to (-,-) R , 
H x H is M-vector space of dimension 4n. On H R the two inner products coincide. Because 
of this we shall simply denote (£,,r]) H or (£,77) the scalar product on H x H whenever 
£,77 G if M . Furthermore, j is norm preserving ||a;|| H = ||j(a:)|| ffR and (x,y) R = (j(a;),j(y)) HR 
for all x,y £ H. The orthogonal complement of i?K in H x H with respect to (•, -) R is given 
by 

(2.22) Hi = iH R = I r-(x + c(x)), ^(x - c(x)) j ,ieif 
The orthogonal projection onto ifjt with respect to the real structure is given by 

(2.23) P R :HxH^H R cHxH, P R {u,v) = Q(« + c(u)),^(y + c(y)) 

The projection onto the orthogonal complement Hj^ is given by 

(2.24) 

p r = 1-Pr , *£(«,«) = f ^(u-c(u)),-(v-c(v))J =i (—(u-c(u)),—(v-c(v)) 
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Note: 

(2.25) P^ (u, v) = -iP{iu, iv) 
Fix a,b G H and define the linear operator on H, 

(2.26) T a , b :H^H, T ajb (x) = (x, a)b 
Associate the following IR-linear operator on H x H 

(2.27) T ay .HxH^H R cHxH, f afi {u) = (u,j(a))J(b) + (u,}(ia))J(ib) 

Note (j(fe), }(ib)) R = and H R is invariant under the action of T ab . Direct computations 
show the following diagram is commutative: 

T a , b : H ^4 H 

(2.28) H Ij 

rp jj a,b jj 

J-a,b ■ Hr > lis. 

Similarly each symmetric operator in Sym(H) gets mapped into a symmetric operator in 
Sym(H^) of double rank. Denote by r this mapping, r : Sym(H) — > Sym(H R ) that makes 
the following diagram commutative: 

H -^ H 

(2.29) j| |j 

f7 ill! U 

If desired, t(T) can be extended to H x if using the IR-linear scalar product (-, -) R on if x H. 
t is IR-linear but not C-linear. In particular: 

(2.30) T=lx,x}eS 1 >° N=> r(T) = [j(x), j(x)] + [j(ix), j(ix)] G 5 2 '°(ii K ) 

(2.31) T=[x,y]G5 1 ' 1 !=► r(T) = [)(*), j(y)] + [j(ix), j(iy)] G 5 2 ' 2 (ii M ) 

In general an operator in S p,q (H) gets mapped into an operator in S 2p,2q (H R ). Note the map 
r preserves scalar products between Sym operator up to a multiplicative constant: 

(2-32) (T,S) B{H) = 2(r(T),r(S)) B{H&) 

Thus r is a monomorphism (i.e. linear injective morphism) between IR-vector spaces Sym(H) 



and Sym(H R ). We shall discuss in more details the linear map r in subsection 3.2 

Note, in the case H = C n , c(z) = (z\, ■ ■ ■ , z n ) T , and (v,w) = v T c(w), then H^ 
{(real(v), imag(v)) , v G C n } = IR 2ri , and the IR-linear map J has the block form 

J-- 



-i 
i 



where i is the identity matrix of size n. 

We return to the frame set F = {/\, • • • , f m }. We let 

(2.33) $ fc = r([/ fc , M) = [p fc , ^] + [.tyfc, J^I 
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denote the image of Fk = [/&, /&] under r, where y2fc = j(/fc) an element of H-r. The last 
notation we introduce here is the following map on H R : 

m 

(2.34) I : H R -> Sym(H R ) , /(£) = £$ fc [£,£]$ fc 

fe=i 

As we will see later J(£) is related to the Fisher information matrix for the measurement 
model (2.3), and it plays a key role in obtaining a necessary and sufficient condition of 
injectivity for the nonlinear map a. Since 

**[£,£]$* = !$*£,$*£! 

/(£) has the following more explicit form: 

in 

(2.35) 1(g) = ^2 K> M , v k = $fc£ = (6, V^rV^ + (f , Jpk) R Vk- 

k=i 

2.2. Existing Results. We revise now existing results on injectivity of the nonlinear map 
a. A subset Z of a topological space (X, r) is said generic if its open interior is dense. In 
the following statements, the term generic refers to Zarisky topology: a set Z C K nxm = 
K™ x • • • x K n is said generic if Z is dense in K nxm and its complement is a finite union 
of zero sets of polynomials in nm variables with coefficients in the field IK (here IK = K or 
K = C). 

Theorem 2.1. In the real case when H = M. n the following are equivalent: 

(1) The nonlinear map a is injective; 

(2) (^ Th.2.8) For any disjoint partition of the frame set F = Fi U F 2 , either Fi spans 
H or F2 spans H . 

(3) ($\,Th.2.4(2)) For any two vectors x,y G H if x 7^ and y ^ then 

m 

52\(x,fk)\ 2 \«y,fk)\ 2 >0i 

(4) (\§\,Th.2.4(3)) There is a positive real constant oq > so that for all x,y G H , 

in 
(2-36) ^Kx,/ fc )| 2 Ky,/ fc )| 2 >a ||x|| 2 ||y|| 2 

fc=i 

(5) ($\,Th.2.4(4)) There is a positive real constant a > so that for all x G H, 

in 

(2.37) R(x) :=Ys\( x 'fk)\ 2 (-Jk)fk > a I 

fe=i 

where the inequality is in the sense of quadratic forms. 

Additionally, the following statements hold true: 

(1) ([4\ Prop2. 5) If a is injective then m > 2n — 1; 

(2) (^Prop.2.5) If m < In — 2 then if cannot be injective; 
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(3) ([A\Cor.2.7(l)) If m = 2n — 1 then a is injective if an only if¥ is full spark; 

(4) (^Cor.2.6) If m > 2n — 1 and ¥ is full spark then the map a is injective; 

(5) f[4] Th.2.2) If m > In — 1 then for a generic frame ¥ the map a is injective. 

In the complex case the following results are known: 

Theorem 2.2. In the complex case when H = C n the following statements hold true: 

(1) (\Q Th.3.3) If m > An — 2 then for a generic frame ¥ the map a is injective; 

(2) ([1A\) For any positive integer n there is a frame of m = 4n — 4 elements so that the 
nonlinear map a is injective; 

(3) f[38j Corollary 4) If ct is injective then 

{2 if n odd and (3 = 3 mod 4 
1 if n odd and (3 = 2 mod 4 
otherwise 

where (3 = (3(ri) denotes the number of l's in the binary expansion of n — 1. 

(4) The following are equivalent: 

(a) The map a is injective; 

(b) f[38j Prop. 2) 

(2.39) ker(A) n (S 1 ' - S 1 ' ) = {0} 

(c) fjllj, Theorem 4) dimS(u) > 2n — 1 for every u G C n \ {0}; 

(d) fjllj, Theorem 4) S{u) = span^iu) 1 - for every u G C ra \ {0}. 

In the last two conditions, S(u) = span^{fkfk u }i<k<m-i where fkfk u ls seen as a ^ n vector 

inM 2n . 

3. New Analysis Results 

This section contains our injectivity results of the nonlinear map a as well as an in-depth 
analysis of the spaces S p,q . 

Theorem 3.1. Let H be a C-vector space of dimension n, with scalar product (, ) and 
conjugation c : H —¥ H. The following are equivalent: 

(1) The nonlinear map a : H —¥ W 71 , (ot(x))k = \(x, fk)\ 2 is injective. 

(2) There is a positive constant a Q > so that for every u,v G H 

m 

(3.40) ^|(F fc ,[ M ,t;])| 2 >a ||Kt;]||?. 

fc=i] 

where Fk = lfk,fk}- Explicitely, this means: 

m 

(3.41) ^2(real((u,f k )(f k ,v))) 2 > a [||u|| 2 ||i;|| 2 - (imag((u, v})) 2 ] . 
fc=i 

(3) For any £ G H R , £^ 0, rank(I{£)) = 2n-l. 
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(4) There is a > so that for all £ G Hr, £ ^ 7 

(3-42) /(0>«o||el| 2 ^ 

where the inequality holds in the sense of quadratic forms in H^, and 

(3-43) P j- i = i--A-lJ^jQ 

is the orthogonal projection in H-g, onto the orthogonal complement of J£ in H^. 

Remark 3.2. The two constants Oq in (2) and (4) can be chosen to be equal, hence the same 
notation. We we see in the next section, this common constant is related to robustness and 
stability of any reconstruction scheme. 

Remark 3.3. The proof of (3) <^ (4) shows that the optimal bound ao is given by 

(3.44) a° pt = min^H R ,u\\=i a n-i(I(0) 

where a n _ 1 (/(^)) denotes the next to the smallest eigenvalue of /(£). 

Remark 3.4. The choice of the nuclear norm and the square in { 3.^0^ is somewhat arbitrary. 
For any p, q > 1 (including infinity), there is a constant a PtQ > so that 

m 

(3-45) EK F ^ v J>l P ^ a ^IIMII" 

fc=i 

An interesting corollary is obtained in the case when a is restricted to a subspace of H. 
Specifically there is sometime the case the underlying signal is real. The canonical description 
of such a condition is to be conjugation invariant. Let H' denote this set: 

(3.46) H' = {x EH ; c(x) = x} 

Note H' is a IR-linear space, but it is not a C-linear space. Since x is restricted to H' 
it follows the equivalence class in H' of x is given by x = {x, —x}. Consequently the 
appropriate quotient space is given by H' = {{x, —x} , x G H R }. Let 7Ti : H x H — > H and 
7r 2 : H x H — > H be the canonical projections onto factors: 7Ti((u, v)) = u and ^((w, v)) = v. 
Then it is immediate to check that H' admits the following equivalent descriptions: 

(3.47) H' = {xeH; tt 2 (j(x)) = 0} = T 1 ({H R n (H x {0})) 

Note in H 1 , (u,v) is always real for all u,v G H' since (u, v) = (c(v),c(u)) = (v,u). Let 
F = {/i, - - - , f m } be the frame set in H. Let 

(3.48) 9k = m(i(fk)), h k = n 2 (i(f k )) 

with 1 < k < m. Note f k = j _1 ((^ fc , h k )) = g k + ih k and g k , h k G H' . Set: 

(3.49) $' fc = \g k , g k \ + {h k , h k \ G S 2 '°(H') C S 2 '°(H) , 1 < k < m 
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Note: 

(3.50) (F k ,lx,xj) = \(f k ,x)\ 2 = (& k x,x} = (& k ,lx,x}) 

(3.51) (F k ,lu,vj) = real((u,f k )(f k ,v)) = (& k u,v) = (& k ,lu,vj) 

for all x,u,v G H' . Thus the linear map A restricted to Sym(H') can be thought of as 
taking inner products with a family of rank-2 nonnegative operators. We have 

Corollary 3.5. Let H be a n- dimensional complex Hilbert space with scalar product (, ) and 
conjugation c : H — >■ H . Let H' = {x G H ; c(x) = x} be the c-invariant set. The following 
are equivalent: 

(1) The restriction to H' of the nonlinear map oc\h> is injective on H' . 

(2) There is a constant a > so that 

m 

(3.52) ^|($^>| 2 >a |M| 2 |M| 2 

fc=i 

(3) For all u^O, 

(3.53) dim R span R {(u,g k )g k + (u, h k )h k ; 1 < k < m} = n 

(4) For any x ^ 0, 



n 



(3.54) ronife[^$'fc[x,x]$'* 

\fc=i 

(5) There is a constant a$ > so that 

in 

(3.55) I'(x) = J2 &klx,x}& k > a ||^|| 2 l^ 

fc=i 

in the sense of quadratic forms on H' , where \n> is the identity map on H' . 



3.1. Analyis of sets S p,q . In addition to (2.11) the sets S p,q introduced in (2.7) have the 
following properties summarized in the following lemma: 

Lemma 3.6. (1) For any p\ < p2 and q\ < q2, S pi,qi C S P2,q2 

(2) For any nonnegative integers p, q the following disjoint decomposition holds true 

v <? 

°r,s 



(3.56) S p ' q = |J |J S 



r=0s=0 



where by convention S 0,0 = S 0,0 = {0}, and <S p ' q = for p + q > n. 
(3) For any nonnegative integers p, q, 



(3.57) - S p ' q = S 
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(4) For any bounded T G B(H) and integers p, q, 

(3.58) T * S P,q T c &>,* 

The inclusion becomes equality if and only if T is invertible. 

(5) For any integers p, q, r, s, 

(3.59) S p ' q + S r ' s = S p ' q - S s ' r = S p+r ' q+s 

Proof of Lemma 13.61 

First three assertions are trivial. Assertions (4) & (5) state that the only invariant of a 
quadratic form over C is its signature. Assertion (4) can be shown using the singular value 
decomposition of T and the continuity of spectrum with respect to matrix entries. For (5), 
one can see this statement as a special instance of the Witt decomposition theorem, a much 
more powerful tool in the theory of quadratic forms especially in the case of vector spaces 
over fields of charactersitics other than 2, see [46]. However for the benefit of those who 
prefer a more direct proof, here is a sketch of this result. Using spectral decomposition and 
rearranging the term, one can easily see that S p+r,q+s C S p ' q +S r,s . For the converse inclusion, 
we need to show that if T, S G Sym(H) are so that T has at most p positive eigenvalues and 
q negative eigenvalues, and S has at most r positive eigenvalues and s negative eigenvalues, 
then T + S has at most p + r positive eigenvalues, and q + s negative eigenvalues. Without 
loss of generality we can assume T G S p,q and S G S r,s . Using spectral decompostions of T 
and S we obtain 

(V p+q \ / r r+s \ 

j2 a kA k - J2 a ^ + [J2 bkBk - J2 hkBk )=u-v , 
k=l k=p+l / \fe=l k=r+l J 

p r q r+s 

u=J2 akAk + J2 hkBk - ° ' v = J2 akAk + Yl bkBk - ° ' 

fe=l fc=l k=p+l k=r+l 

where A k , B k are rank one orthogonal (spectral) projectors with A k Aj = and BiBh = for 
all k 7^ j and I ^ h, and a k , b k > 0. Thus U G S p+r, ° and V G S q+S, °. The claim now follows 
by induction provided we show that for any R G S a,b and E G S 1,0 a rank one associated to 
vector e G H (E = (■, e)e), then E + Re S a+1 > b . 

Indeed this last asserition is shown as follows. Let R = J2t=i c k('i 9k)9k — Ylk=i+a c k('i 9k) 9k 
be its spectral decomposition. Two cases are treated distinctly. 

Case 1: If e ^ R(H) then P* = {^i, • • • , g a +b, e} is linearly independent. Let T = {7^ , < 
k < a + b + 1} be the (unique) biorthogonal system to V. Let T be an invertible operator 
that maps an orthonormal set A = {5i, ■ ■ ■ , S a+ b + i} into T, T5 k = 7fc, l<k<a + b+l, 
and maps the orthogonal space to A onto the orthogonal complement to T. Note T*g k = 5 k , 
1 < k < a + b and T*e = b~ a +b+i- Then a direct computation shows that 



T*(R + E)T = Y, c k (; 5 k )5 k + (•, S a+b+l )S a , 



2+b+l 
k=l 
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Thus the spectrum of T*(R + E)T is composed of {ci, ■ ■ ■ ,c a , — c a +i, ■ ■ ■ , c a+ b, 1} which shows 
that T*{R + E)T G S a+1 ' b , and by (4), R + E G S a+1 > b . 

Case 2: e G R(H). The rank of R+E is the same as the rank of R, hence dim ker(R+E) = 
dim ker(R). Hence R + E G S a ,b with a' + b' = a + b. Now by the continuity of spectrum 
with respect to small perturbations, it follows that, for a small perturbation e h- > e' $. R(H), 
R + E' G S a "' b " with a" > a', b" > b'. But the proof of case 1 shows S a "' b " C S a+1 ' b . Hence 
b > b" > b' > b, and a + 1 > a" > a' > a. Thus R + E G 5 a+1 ' 6 . Q.E.D. □. 



Next we analyze the space S 1 ' 1 in more detail. The special factorization we obtain here 
can be extended to other classes S p,g but we do not plan to do so here. First set the following 
matrices: 



(3.60) 



H 



1 

1 



D 



1 








V 



1 
7! 



i i 

-i i 



In the next result we use a generalized unitary group U(l, 1; H). Recall its definition. 
Definition The groups U(l, 1) and U(l, 1; H) are defined by 

(3.61) U(l,l) = {AeC 2x2 , A*DA = D} 



(3.62) 



U(l, 1; #) = {A G C ix2 , A*iM = H} 



These groups been studied extensively in literature. See for instance [56], section 10.4. In 
particular the two groups above are unitarily equivalent to each other, and the matrix V 
provides such an equivalence: 



(3.63) 



A G £7(1, 1;H)<*B = VAV* G U(l, 1) 



The quadratic form u)(zi, z 2 ) = \zi\ 2 — \z 2 \ 2 is invariant under the action of U(l, 1), whereas 
4>{z\, z-i) = Z\Z 2 + z\Zi is invariant under the action of U(l, 1; H). 

Lemma 3.7. (1) S 1 ' 1 = S lfl - S lfl = 5 1 ' + S ' 1 . 
(2) For any T G S ,x there are u,v G H so that 

1 



(3.64) 



T 



((-,u)v+ (-,v)u) = {u,vj 



(3.65) 



where the symmetric outer product \u,v\ has been defined in (2.5). 

If T = Oi[ei,ei] — a2[e2,e2] with 01,02 > and (ek,ej) = Skj is its spectral 
factorization then 

u = 



Oiei 



fl2e 2 , v 



a,\e\ 



a 2 e 2 



provides a particular factorization in (3.64). Conversely, for any u,v G H, the 
operator T = [w, vj is in S 1 ' 1 . 
(3) Let T = [ttojVo] G "S 1 ' 1 . Then any pair (u,v) of vectors in H so that T = {u,vj is 
given by 

(3.66) 



u 



anUQ + a 12 v , v = a 21 u + a 22 v 
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for some matrix A = (ak,i)i<k,i<2 with A e [7(1,1; if). Conversely, for any matrix 
A e U(1,1;H), {u,vj = luo,voJ where (u,v) are given by (3.66). 
(4) Let T = [w, vj. Then traces and spectrum Sp(T) = {a + ,a-} are given by 



(3.67) 
(3.68) 

(3.69) 
(3.70) 

(3.71) 



tr{T} 
tr{T 2 } 



a + 



real((u, v)) = (u, v) 

1 

4 

1 

2 

1 
2 



(((u,v)) 2 + ((v,u)) 2 + 2\\u\\ 2 \\v\\ 2 

2 II II 2 i / \2 /• \ 2 \ 

H +( u , v )r- (*m,«)kJ 
(W)«)m+ vlMI 2 NI 2 + (*«>«)» 



( m ^)r- yll 

Hence the nuclear norm of T is given by 

nil 



n2 n n2 



(z«,u)i 



>o 

<0 



a i + a_ 



/ll l|2|i i|2 i /■ \z 

y ||u|| ||v|| + {iu,v) I 



Remark 3.8. 1. The need for studying S 1,1 arose from the bahavior of the IRLS algorithm 
described in [H] . However it quickly become apparent the space S 1,1 and its factorization given 
by (3.64) are crucial for understanding the injectivity of the nonlinear map a especially in 
the light of Theorem \k~2] (3), equation (2. 3 Sty . 

2. The choice in (3.65) has the following two additional properties: 



(3.72) 
(3.73) 



ll u o|| 
(u ,vq) 



a± — a 2 



\T\ 



\/ai + 02 = V a + ~ a - - 

= a + + a_ = tr{T} (a real number!) 



where HT^ represents the nuclear norm of T (sum of its singular eigenvalues), and a\ 
a + > and 02 = —a- > are its singular eigenvalues. 

Proof of Lemma 13.71 



(1) is a direct application of Lemma 3.6 [5). 

(2) Since [, ] is IR-linear and [ei,e2] = [e2,ei] we obtain 

lu ,v j = ai[ei,ei] + y/a 1 a 2 le 1 ,e 2 j - y/a 1 a 2 le 1 ,e 1 j -a 2 [e2,e 2 ] = ai[ei,ei] -a 2 [e2,e 2 ] = T 

For the converse, note rank(T) = dimspan{u,v} < 2. Then the claim follows once we 
establish (3.69) and ( 3.70[ ) which we do below. 

(3) A direct computation shows 

K^l = ^u^i + a 21 a n )(-,u )u + -{aua 22 + a 21 a 12 )(-,u )v 

1 1 

+ ^,(a 12 a 21 + a 22 a 11 )(-,v )u + -(ai 2 a 2 2 + 022012) (-,^0)^0 
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Since {uo,vo} are linearly independent, [it, vj = |[ito,i>o] implies 

a n a 2 i + a 2 iOn = , ana 2 2 + a 2 ia 12 = 1 

^120-21 + ^22^11 = 1 , ai 2 a 2 2 + ^22^12 = 

Which corresponds to A* HA = H. Hence A G 17(1,1; H). Conversely, if A G U(1,1;H) 
then above relations are satisfied which imply [u, vj = [wo,fo]. 

(4) The equation (3.67) comes from the definition (2.5) and the fact that tr((-,u)v) = 
(v,u). For T 2 compute first 

T 2 — j((v,u}(-,u)v+\\u\\ (■,v)v + \\v\\ (•,u)u+(u,v)(-,v)u) . 



Then (3.68) follows from this relation and real(({u,v}) ) = (real((u,v))) — (imag((u,v))y 
Finally, (3.69) and (3.70) come from solving: 



(3-74) ., „, ,,,,.„ 



a + + a_ = tr(T) 
a 2 + + a 2 _ = tr{T 2 



Q.E.D. a 



Topologically, S 1,0 and S 1,1 are not differentiable manifolds. Instead they are algebraic 
varieties, since are given by the zero loci of certain polynomials. We have the following 
result: 

Lemma 3.9. (1) The set S 1,0 is an analytic manifold in B(H) of real dimension 2n — 1. 
As real manifold, its tangent space at X = [x,x] is given by 

(3.75) T x S lfi = {lx,yj , y e H} 

The M.-linear embedding H \-¥ TxS 1 ' given by y h-> (p x (y) = {x, y\ has null space given by 
ker (p x = {iax ; a G M.}. 

(2) The set S 1,1 is an analytic manifold of real dimension An — 3. As real manifold, its 
tangent space at T = {x,yj is given by 

(3.76) TVS 1 ' 1 = {lx,uj + ly,vj , u,v G H} 

The M-linear embedding H x H i-> TVS 1 ' 1 given by (x,y) h- > (p Xjy (u,v) = {x,uj + {y,v} has 
null space given by ker (p x>y = {a(ix, 0) + 6(0, iy) + c(y, —x) , a, b, c G K}. 

Proof of Lemma 13.91 

Let Ci, • • • , c n : Sym(H) — > R be the coefficients of the characteristic polynomial: 

de£(s/ - T) = s n + c^s™" 1 + c 2 (T)s n - 2 + • • • + c n (T) 

Thus Ci(T) = -tr(T) and ^(T) = (-l)"c/et(T). 

(1) The manifold structure can be shown as follows. First note that 

S + = {S e Sym(H) ; c x (T) = -tr(S) < 0} 
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is an open subset of Sym(H) and therefore a manifold of same real dimension as Sym(H) 
(which is 2n). Next note 

Hence S 1,0 is an algebraic variety. We obtain further it is an analytic manifold by showing 
its tangent space has same dimension at all points. Let X = [x ,ico] G S 1 ' . We consider 
the set of all differentiable curves 7 : I — > <S 1,0 , 7(0) = X, G I C M. open, passing through 
X. Then the tangent space at X is given by 

TxS 1 ' = {j tl (t)\t=o} 

For each such curve, 7 : I — > (S) 1 ' there is a unique curve x : / — >■ H such that j(t) = 
lx(t),x(t)J with x(0) = Xq. In fact, locally, 

x(t) = 7(t)(x ) 

which shows x(t) is differentiable. Then a direct computation shows 

j t l(t)\t=o = [x(0),x(0)] + [x(0),i(0)] = 2[x ,x(0)] 

Sin ce i;(0 ) can be chosen arbitrarily in H, it follows the tangent space to S 1,0 at X is given 
by (3.75). The real dimension of the M-vector space TxS 1,0 is 2n — 1 once we notice the 



kernel of the IR-linear map ip x is one dimensional and given by the span of ix. 
(2) The algebraic manifold structure is given by the intersection 

5 1 - 1 = c 3 - 1 (0)n---nc; 1 (0)n5- 

where 

S— = {S e Sym(H) ; c 2 {T) < 0} 

The analyti c man ifod structure follows once we establish the tangent space at every point 



is given by (3.76) and has same dimension An — 3. Fix T e S 1,1 , T = [xo, 2/0] with {xo,yo} 
linearly independent, and let {7 : J — > S 1 ' , G I C R open , 7(0) = T} be the set 



of differentiable curves passing through T. By Lemma pi we know 7(2) = [x(t),y(t)J for 



some x : I — ^ H and y : I — > H . However these functions are not unique. However we can 



choose them to be given by the spectral factorization of 7(f) via (3.65). Furthermore a direct 
application of functional calculus (see [53]) shows that T — > P\ G 1S 1 ' and T ^ P 2 £ S 1,0 
are analytic, where T = P\ — P 2 is its spectral decomposition with P\P 2 = 0. Hence 
£ G J H- -Pi(i) and t £ I ^ P 2 (t) are differentiable. Then the component functions x(t) and 
y(t) now uniquely defined by: 

(3 - 77) x{t) = Tfwmrm m Pimxm ' y{t) = Ttwmrm m P2mm 
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are different iable. The derivative at t = is given by 

j t l(t)\t=o = 2[x(0), |/(0)] +2[y(0), i(0)] 

And since (x(0),y(0)) and (x ,y ) are related by a [7(1,1; if) matrix, and the fact that 
(x(0),?/(0)) is arbitrary in if x if , we obtain the tangent space is given by (3.76). A direct 
computation shows that (ix, 0), (0, iy) and (y, —x) are independent vectors in the null space of 
<f X; y which concludes the proof. □ 

Remark 3.10. Compared to the complex projective manifold CP n , S 1,0 is diffeomorphic to 
CP n x ]R + . The extra IR + component comes from the fact that rank-1 operators in S 1,0 have 
arbitrary trace. This explains the real dimension of S 1,0 : 2dini]KCP n + 1 = 1(n — 1) + 1 = 
2n — 1. On the other hand, using spectral factorization, S 1,1 has real dimension given by 
dim K 5 1 ' + dim K 5 1 ' — 1 = 2(2n— 1) — 1 = An — 3. The —1 term comes from the orthogonality 
of the two eigenvectors of an operator in S 1,1 . 



Remark 3.11. As suggested by Bernhard Bodmann [To], it can be shown the subset of 
projections inside S 1,0 is in fact a Kahler manifold (diffeomorphic to CP n ). A similar 
conclusion holds for S 1 ' 1 since the set of projections inside S 1 ' 1 is the same to the set of 
projections inside S 1,0 . However none of these two manifolds, S 1,0 and S 1 ' 1 , is Kahler: they 
have odd real dimension, hence cannot admit a symplectic structure. 

3.2. Analysis of the linear map r . We introduced earlier the M-linear map r that maps 
Sym(H) operators into Sym(H R ) operators, using the real linear structure on these spaces. 
In order for the diagram (2.29) to be commutative, the map r(T) is given explicitely by 

(3.78) r(T)(£)=j(T(j- 1 (0), f G H R 
The followin lemma summarizes the basic properties of the map r. 

Lemma 3.12. (1) Let P be an orthogonal projection of rank k in H. Then t(P) is an 
orthogonal projection of rank 2k in if jr. Furthermore if {ei, ■ ■ ■ , e^} is an orthonormal 
basis in the range of P, then {j(ei), - • • , j(efc), Jj(ei), ■ • - , Jj(efc)} is an orthonormal 
basis in the range ofr(P). 

(2) If T G Sym(H) has spectrum (ai, a 2 , ■ • • , a n ) then t(T) in B(H R ) has spectrum 
(a l5 ai, 02, a2, ■ ■ • , a n , a n ). 

(3) For any two operators T,S G Sym(H), r(T),r(S) G Sym(H R ) and 

(3.79) tr{r(T)r(S)} = (r(T), r(5)> fl(flk) = 2(T, S) B{H) = 2tr{TS} 

(4) Let 1 < p < oo. The p-norms of a symmetric operator T G Sym(H) and t(T) G 
Sym(H^), are related by 

(3-80) ||r(T)|| p = 2 1 /'||T|| p , if p < oo 

(3.81) ||T|| = ||t(T)|L = ||T|L=||T|| 
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Proof of Lemma 13.121 

(1) First we prove the statement for rank-1 projectors. This comes from directly checking 
equation (2.30). Thus P = [e, e] gets mapped into r(P) = [£, £] + [J£, J£], where £ = j(e). 
Next if {ei, • • • , e k } is an orthonormal basis for the range of P then 

k 

P = J2lei,e l j 
i=i 



By IR-linearity, t(P) has the form 

k 

i=i 

where £/ = j(ej), 1 < I < k. Note 

(J£i,J£ s ) = (6,6) = (ei,e s ) M = <5 M , (J6,6) = (*ej,e 8 ) B = 

Hence {£i, •••,£&, J£i, • • • , ■/£&} is an orthonormal set and since it is spanning the range of 
t(P) it follows it is an orthonormal basis in the range of r(P). Hence r(P) is an orthonormal 
projection on H-^ of rank 2k. 

(2) Follows by using the spectral factorization of T, 

s s 

T = J2 a kPk H- r{T) = Y, a k r{P k ) 
i=i i=i 

Then r(P k )r(Pi) = for k ^ I. Hence the equation for t(T) above represents its spectral 
factorization, the conclusion follows. 

(3), (4) Follow from (2). □ 



3.3. Proof of Theorem 3.1 and its Corollary 3.5 (1)<^>(2). According to Theorem 
|2.2[ 4), nonlinear map a is injective iff ker(A) fl (S 1 ' — S 1 ' ) = {0}. But using Lemma 3.7[ l) 
and (2) we get equivalently that a is injective iff for all u,v G H with [u, vj ^ 0, 

-4(M)^o 

Equivalently this means 

m 

Y,\(F k ,lu,v})\ 2 >0 



k=l 



fl,l 



Consider now the unit ball in S 1 ' 1 with respect to the nuclear norm, say S^' 1 . This set is 
compact in Sym(H). Then let 



(3.82) 



a = min TeS i,i m ^ =l ^ | (F k , T}\' 



k=i 



By homogeneity we obtain (3.40). Then 



(F k , {u,v}) = -((uj k }(fk,v) + (v,f k )(f k ,u)) = real((uj k )(f k ,v}) 
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and by Lemma |3.7[4) 



fr tiii2 ii ii2m i 

\\u,v INI, = \\u\\ V 



(iu,v) R = \\u\\ \\v\\ — (imag({u,v))) 2 

Putting together all previous derivations we obtain the equivalence (1) -v4> (2). 
(2)<^(4). Using Lemma 3.9[3) we obtain (3.40) is equivalent to 



(3.83) 



J2 \(r(F k ),r(lu,vj))\ 2 > 4a [\\u\\ 2 \\vf - (real((iu,v))?] 



fc=i 



Now ( |2~30| and fl2~3l| ) imply 

r(F k ) = \ip k , ip k l + \Jip k , Jip k j = y k , ip k \ + Jfipk, ipkjJ* 

r(lu, v}) = U, rij + [J£, Jrj\ = It, rjj + Jit, V p* 
where ip = j(f k ) and £ = j(w), rj = j(i> ). A direct computation using J* = —J shows that 

(T(F k ) 1 lJt,J V }) = (T{F k ),lt,Tl}) 

Thus 
(r(F fc ),T([«,u])) = 2(l<p k ,<p k } + lJ<Pk,J<Pkl,U'Vl) = 2 [(^,fk)(fk,v) + (£,J<Pk)(J¥k,v}] 



With the equation (2.35) we obtain 

m 

(r(F fc ),r([«,T;])> = 2(v klV ) =»£|<t(F*),t(|[- 



w, f 



4</(0»?,»?> 



fc=i 



Now the right-hand side of (3.83) is processed as follows. Note ||w|| = ||£||, \\v\\ = \\r]\\, and 



Thus 



real((iu,v)) = (iu,v) M = (J£,r)) 



\u\\ 2 \\vf - (real(( lU ,v))) 2 = Ufllvf ~ ((^» 2 = ((ll^fl " 1^, J®)**) 



Substituting in (3.83) we obtain (3.42). 



(3)-t=r-(4). Assume rank (1(0) = 2n — 1 for all £ ^ 0. A direct computation shows that 
I(0(J0 = 0. Hence J£ is the only independent vector in ker(I(0)- It follows there is 
a = a(£) > so that 

1(0 > a(0PJs 
Note the a(£) represent the smallest nonzero eigenvalue of I(£) which must be the 2n — 1 th . 
Using the homogeneity of degree 2 on the unit sphere of Hr, we get a(£) = ||£|| oCnfir) > 

a ||£|| which proves (3.42). Conversely, if ( J3.42 ) holds true, then J(£) has rank at least 2n— 1. 
Again since J£ is in the kernel of /(£), it follows that /(£) must be of rank exactly 2n — 1. 
Q.E.D. D 



Proof of Corollary 3.5 
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(1)4^(2) follows from Theorem |3.1 [ 2) and equation (3.51 ) and the fact that imag((u, v)) = 
for all 11,1)6 H' . 

(2)<^>(5) follows from the relation 

m m 

fe=l fc=l 

and 

m m 

fc=i *;=i 

(4)4^(5) follows by the homogeneity of I'(x) with resepct to ||x||. 

(2)<^>(3) follows by contradiction and homogeneity of (u, v) \-± Y^k=i \(&'kU,v)\ 2 with 
resepct to the product ||u||||t?||. 

Q.E.D. D 

4. Performance Bounds on Reconstruction Algorithms 

In this section we present two performance bounds applicable to any reconstruction algo- 
rithm. One bound is deterministic and is based on the constants a® introduced in Theorem 



3.1 The second bound represents the Cramer- Rao Lower Bound for the stochastic model 



(2.3) 



4.1. Lipschitz bounds of the inverse map. Consider the nonlinear map a : H — > IR m . 
We shall establish a deterministic performance bound for any inversion algorithm in terms 
of the Lipschitz bounds of the map: 

(4.84) A : S 1 ' ->• E m , A{{x, xj) = a{x) 

Specifically we want want to bound above and below the following expression: 

(4.85) R{Xiy)= \\M*,4)-Ah,«m\ 

||[a;,a;] - [[2/ , 2/1 1 1 ! 
Since [x, x] — {y, y\ = {u, v} G S 1,1 for some u, v G H it follows: 

(4.86) ,K«] = sup Sk^k^ 

(4.87) inf R(x ,y } = inf ^^1 

x,yeH u,veH ||[m, ^Jllx 

These ratios can be further processed as follows 

£r=il(^M>l 2 __ (i(Z)v,v) 



Mil? " U\\ 2 (p^v) 
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where £ = j(w) and 77 = j(i>). Since I(£)?7 = HO-PJ^V ^ follows: 

JA0^L = M^ = max || /(0 n 

and 

inf (i(Ov,v) = a o P t 
t*#>\\Z\\ 2 {Pj?i,v) " 

Thus we proved: 

Theorem 4.1. Assume the nonlinear map a is injective. Then the map A : <S 1,0 — > 



defined in (4-84) is bi-Lipschitz between (iS 1,0 , || • HJ and (R m , || • ||) m£/i the Euclidian norm, 



and it has the upper Lipschitz bound 
(4.88) 

and the lower Lipschitz bound 



(4.88) B = max >/P(0F 

£e-ffK,||£||=i 



(4.89) A„=^r = ^ 6 min M0 „_ l( /K)) 

Specifically 

(4.90) A)IIM - MIL < M(M) - *4(M)|| < B |||[s,ar] - [y,y]| 
/or a// u,v £ H. 



4.2. The Cramer-Rao Lower Bound (CRLB) . Consider now the noise model (2.3[) 



We are interested in obtaining a lower bound for any unbiased estimator of x. The derivation 
of the CRLB in this paper coincides with the one presented in [TT] (see Theorem 23 there). 
In turn this follows the receipe presented in [D]. We will just present the key steps of this 
derivation. Note that our derivation is canonical, that is basis independent. 

Due to non-holomorphy of the nonlinear map a, the analysis is done in the realification 
space Hjh. We denote £ = j(x) for the signal x 6 H. The frame set if F = {/1, • • • , f m } 
and Fk = [/&,/&] G S 1,0 (H) denote the measurement operators. Recall our notation $& = 
r(Ffc) = [v?fc,(y9fc] + {Jtpk, JipkJ £ <S 2 '°(Hr). First we compute the Fisher information matrix 
associated to £. The likelihood for this problem is 
(4.91) 

1 / 1 \ 1 ( I m 

p(ylO = (2^«? [-^\\y - «(*)II 2 J = (27r)WVm ^ (^"2^ l2/fc " ($fee ' e)l 

The Fisher information matrix is given by (see [45] ) 

(4.92) F(£) = E [(V c log(p( 2 /|0))(V ? log(p(y|£))) T ] 
The canonical form of this operator is 

(4.93) F(f) = E[[V e log(p(y|0), V e log(p(y|0)]] 
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A little bit of algebra shows 

A rn A 

(4.94) F(0 = — J2* k lt,tl*k = V(0 

fc=i 

In general the covariance of any unbiased estimator is bounded below by the inverse of the 
Fisher information matrix (operator). However in this case the Fisher information operator 
is not invertible. This fact simply expresses the statement that x is not identifiable from the 
measurements y G M. m alone. As we know the nonlinear map a is not injective on H but 
instead it is injective on H . The nonuniqueness on H is reflected in having a singular Fisher 
information matrix on if jr. To solve this issue we need to fix the global phase factor. One 
solution is to fix a basis and decide that a particular component (say the last component) 
is real. Such an approach was taken in [11]. Here we propose a canonical solution to this 
normalization. An oracle provides us with a vector z G H, so that (x, Zq) > is positive 
real. Assume Zq is normalized ||zo|| = 1- Note there are two pieces of information that can 
be extracted from here: First the fact that x is not orthogonal to z ; in particular x ^ 0. 
Second, the global phase to recover x from its rank-1 operator [x, x\ is uniquely determined 
by the fact that imag((x, Zq)) — and real({x, Zq)) > 0. 

Under this scenario we want to analyze the Fisher information operator obtained earlier. 
Let ipo = j(z ) G H R . We know 

(f , V>o) = real({x, z }) > , (f , Jipo) = imag({x, ip )) = 

Let n denote the orthogonal projection onto the complement of Ji/jo, 

(4.95) n : H R -> H R , U = 1 - [J^o, ^o] 

and let E C H R denote its range, E = {J-^o}" 1 - The oracle tells us £ G E. Hence the class 
of estimators for £ should include only functions 

(4.96) u : R m ->■ j _1 (-B) = {x G H ; (x, z ) > 0}) 
In this case the appropriate Fisher information operator should be 

A rn 

(4.97) f(0 = nF(on = - Vn$ t [u]$ t n 

fc=i 

The following lemma proves that under the scenario described here F(£) is invertible on 
E. 

Lemma 4.2. Assume a is injective on H and zq G H is so that (x,zq) > 0. Then 

(4.98) f(0 = iLF(0n > 4 a ol(^^o)l 2 n 



(7 



where a is the same lower bound introduced in Theorem \3.1\ whose optimal value is given by 

Km. 
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Proof 



We know I(£) > ao||£|| Pj^ from Theorem 3.1 (4). We only need to show 

||£|| 2 np^n> |(^o>l 2 n 

and the inequality is tight. Without loss of generality we can assume ||£|| = since all expres- 
sions are homogeneous in ||£||. Then we need to show that for any u G E, \\u\\ = 1, 

(Pjsu,u) >|<£,Vo>| 2 
This follows from 



inf l-|(u,J£)| 2 = l- max |(u,J£)| 2 = 1 

||u||=l,ueS \\u\\=l,u€E 



(^- jo 2 



k^p 



Note also the inequality is tight. Q.E.D. □ 

Thus we established that F(£) is invertible on E. See also [TT], Lemma 22, for a similar 
statement. Let H ZQ denote the following open set; 

(4.99) H Z0 = {xeH , (x,z )>0} 
Recall an estimator u : M m — > H is said unbiased if 

(4.100) E[w(y)|x] = x 
We can now state the main result of this section: 

Theorem 4.3. Assume the nonlinear map a is injective and fix a vector zq G H . For any 
vector x G H Zo the covariance operator of any unbiased estimator u : lR m —¥ H Zo of x is 
bounded below by the Cramer-Rao Lower Bound (CRLB) given by 

2 ( m ^ 

(4-101) Cov[u(y)] > °- \Y,*Ai,®Pi 

\k=i 

where t denotes the pseudoinverse operator. In particular the Mean Square Error of u is 
bounded below by 

2 [ / m 

(4.102) MSE(co)=E[\\x-u(y)\\ 2 \x}>^trl (£>*[£,£]$ 

The previous Lemma allows us to to predict an upper bound for the MSE of any efficient 
estimator (that is an unbiased estimator that achieves CRLB): 

Corollary 4.4. Assume u : W m — » H Zo is an unbiased estimator that achieves the CRLB 
(4-101). Then its Mean-Square Error is bounded above by 

2| n ^ (2n-l)a 2 




(4.103) MSE(u)=E[\\x-uj(y)\\ z \x]< 



4ao Pi |(x,z )| 5 
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5. The Regularized Iterative Least-Square Algorithm 



Consider the additive noise model in (2.3). Our data is the vector y G M. m . Our goal is to 



find aiGif that minimizes \\y — a(x)\\, where we use the Euclidian norm. Set 

m 

(5.104) J (X) = J2 \V* - (Xfk,fk)\ 2 , Jo ■■ Sym(H) -)> R 



fc=i 



and notice J ([x,x]) = \\y — cx(x)\\ . The least-square error minimizer represents the Maxi- 
mum Likelihood Estimator (MLE) when the noise is Gaussian. In this section we discuss an 
optimization algorithm for this criterion. 

Consider now Jo = \\y — A(X) || restricted to S ,0 , which is an analytic manifold. Consider 
a current point X* = [x*, x*] in an iterative process. Then a descent direction can be thought 



of as a vector in the tangent space to the manifold. According to Lemma 3.12 (1), the tangent 
space at X* is given by operators of the form [a;*, 5J. Since X 1 + [x*, 5] = [x*, x* + 5} G S 1,1 
, one would need to project X 1 + [x*, S\ back into <S 1,0 and choose direction S that minimizes 
(or at least decreases) J (P(X t + [x 4 , 5})), where P is the (nonlinear) projection in Sym(H) 
onto iS 1 ' . However since J is well defined on S 1 ' 1 we choose to optimize S without projecting 
back onto S 1 ' . Thus we obtain the iterative process: 

x t+1 = argmiriuJodx 1 ,^) 

However this process is not robust to noise, the main reason being ill-conditioning and 
multiple local minima of Jo o n S 1 ' . Instead we choose to regularize this process and thus 
to introduce a different optimization criterion. 
Consider the following functional 

(5.105) J : H x H x R + x R+ -> M + 

m 1 

J(u,v,\,n) = ^2\yk- 7j((%fk){fk,v) + (v, fk)(fk,u))\ 2 + X\\u\\ 2 + /j,\\u ~ v\\ 2 + X\\v\\ 2 . 



k=i 



2 



Our ultimate goal is to minimize Jq{u) = \\y — a{u)\\ = J(u,u, 0,/x) over u, for some (and 
hence any) value of \x G 1R + . Our strategy is based on the following iterative process: 

Step 0. we initialize x° as the global optimal solution for a specific pair (A ,/A))- 

Step 1. Repeat: 

(5.106) x' +1 = argmiriuJiu.x 1 ,\ tl Ht) 

Update A 4+ i,/i t+ i according to a specific policy; 
Step 3. Stop when some tolerance level is achieved. 



As we describe below the update (5.106) can be tweeked o achieve a more robust behavior 



(see ( |5l24| ). 

5.1. Initialization. Consider the regularized least-square problem: 

min u J(u, u, X, 0) = min u \\y — oi(u)\\ + 2A||m|| / 
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Note the following relation 



in 



J(u,u,X,0) = \\y\\ 2 + 2X\\u\\ 2 -2j2yk\(uJk)\ 2 + J2\^f^ 1 

fe=i fc=i 

m 

(5.107) = \\y\\ 2 + 2((\I-Q)u,u) + J2\(u,fk)\ 4 



where 

m m 

(5.108) q = J2y k ifk,fk] = J>* F * 



fc=l fc=l 



For A > ||<5|| the optimal solution is u = 0. Note that if Q < (as a quadratic form) 
then the optimal solution of min u \\y — <p(u)\\ is u = 0. Consequently we assume the largest 
eigenvalue of Q is positive. As A decreases the optimizer remains small. Hence we can neglect 
the forth order term in u in the expansion above and obtain: 

J{u, u, A, 0) « \\yf + 2({XI - Q)u, u) 

Thus the critical value of A for which we may get a nonzero solution is A = maxeig(Q) 
is the maximum eigenvalue of Q. Let us denote by Oi this (positive) eigenvalue and e\ its 
associated normalized eigenvector. This suggests t o initia lize A = pa x for some < p < 1 
and x° = (3ei, for some scalar (3. Substituting into (5.107) we obtain 

in 

J(/3e 1 ,/3e 1 ,pa 1 ,0) = || 2 /|| 2 -2(l-p)a 1 /3 2 + (^|(e 1 ,/ fc )| 4 )/3 4 



fc=i 



For fixed p, the minimum over /3 is achieved at 



(5.109) fa = J i 1 |/ )a !,, 4 , x° = /3ei 

V l^k=i K e i)/fe/l 

The parameter /i controls the step size at each iteration. The larger the value the smaller the 
step. On the other hand, a small value of this parameter may produce an unstable behavior 
of the iterates. In our implementation we use the same initial value for both A and p: 

(5.110) yU = A =pai 



5.2. Iterations. Minimization (5.106) is performed in the space H-&. Let £' = j(x*) and 
V = J( M )- Then 

m 

(5.1H) j( u , x\ A*, p t ) = J2 \(i^e, $k?h, v) - y k \ 2 + a,i|77h 2 + p t \\ v - e\f + x t \\ew 2 

fc=l 
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Note the criterion is quadratic in rj. The unique minimum is given by solving the linear 
equation: 

(m \ / m \ 

for £' +1 . In our implementations we decrease (A t ,/i t ) but we limit fi t to a minimum value. 
Thus our adaptation policy is 

(5.113) A m = 7 At 

(5.114) fjLt+i = max(7/i t ,/i mm ) 
where < 7 < 1 is the rate parameter. 

5.3. Stopping Criterion. In our implementation we repeat the iterations until A reaches 
a preset value \ mm . As proved later in this section, the error is linearly dependent on A. 

Alternatively, one can stop the iterations once the modeling error becomes of comparable 
to noise variance. Specifically, a stopping criterion would be 

m 

(5.115) ^b fe -|(^,/ fe )|T<^^ 2 

fc=i 

where k > 1, for instance k = 3. 

5.4. Convergence and Optimality. Consider the following three functionals Ji, J2, J3 : 

Sym(H) xl + xl + 4l defined by 

m 

(5.116) Ji(X,A,/i) = ^|y fc -(X,F fc )| 2 + 2(A + ;U )||X|| 1 -2/itr{X} 

fe=i 

(5.117) J 2 (X,A,/i) = ^|y fc -(X,F fc )| 2 + 2Aa ma:c (X)-4/io min (X) 

fe=i 

m 

(5.118) J 3 (^,A,/i) = ^|^-(X,F fc )| 2 + 2A||X|| 1 -2(2A-/i)a mm (X) 

fe=i 

where a max (X) and a m j n (X) are the maximum and minimum eigenvalue of X, respectively. 
We can prove the following result: 

Lemma 5.1. (1) When restricted to S 1,1 the three criteria coincide: 

(5.119) J 1 (X,\, fJ ,) = J 2 (X,\,n) = J 3 (X,A,/i) , VlG5 1 ' 1 ,A 1 /iGR. 

(2) On Sym(H), the three criteria Ji, J2, J3 are convex. 

(3) T/ie minimum value of J±, J2, J3 on S 1,1 coincides with the minimum value of J on 
H x H: 

(5.120) min J 1 {X,\,n)= min J 2 (X,A,/i)= min J 3 (X,A,/i)= min J(u, v, A, /j,) 

xes 1 ' 1 xes 1 ' 1 xe^ 1 - 1 u,veH 
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for any A, \x > 0. Any minimizer X G 5 1 ' 1 /or J l5 J 2 ? <^3 an ^ (w, "0) for J satisfy 

(5.121) X = {u,vj , \\u\\ = \\v\\ , imag((u,v)) = 

(4) Restricted to S 1,0 all four criteria coincide: 
(5.122) 
J(u,u,X,fi) = Ji([u,w],A,aO = J 2 ([m,w],A,//) = J 3 ([tt,tt],A,/x) = ||y - a(u)\\ 2 + 2A||w|| 2 

and are independent of /i. 

Proof 

(1) and (2) can be checked directly. 

(3) and (4) follow from the following relation: 

J(u,v,X,fi) - Ji({u, v}, A, n) = (A + //)[( ||m || - |H|) 2 + 2||w||||u|| - 2y ||m|| ||u|| 2 - (imag((u,v)))' 
(5.123) > (X + mu\\-\\v\\) 2 



that follows from (3.67) and (3.71 ). 

Remark 5.2. The criterion J 2 shows that the two regularization terms A(||w|| + ||v|| ) and 
fi\\u — v\\ have different effects on the optimizer: the larger the parameter \i the closer the 
lower eigenvalue is to zero; hence the closer the optimizer X is to a rank-1 operator; on the 
other hand, the larget the parameter X the larger the cost of the S 1,0 component in X ; hence 
\\X\\ remains bounded. 

Remark 5.3. In the optimization problem inf Uit)G # J(u,v, A,/i) the constraint is convex but 
the criterion is not jointly convex in (u,v). It is however convex in each individual variable 
u and v. On the other hand in the optimization problems infxesM J S {X,X, fi), 1 < s < 3, 
the criterion is convex in X , but the underlying constraint X e S 1 ' 1 does not define a convex 
set. 

The optimization procedure outlined in the previous subsection describes the following 
mechanism. Consider a path (x*)t>o in H. Set 

X t+1 = \x\x t+1 \ = i ((;X t+1 )x t + (;X t ) X t+1 ) 

Thus a trajectory (x*) 4 >o in H is mappped into a trajectory (X*) t > in S 1 ' 1 . Then the 
algorithm chooses X t+1 along a tangent direction to <S ' at [a;*, a; ], namely a direction of 
the form [x*,w] for some u G H that is of maximum descent for J s , 1 < s < 3. Since the 
algorithm performance is completely characterized by the sequence {X t )t>o, and since for a 
given X G S 1 ' 1 the minimum of J(u,v, A,//) over u,v G H subject to {u,vj = X is achieved 
at a pair (u, v) so that ||u|| = ||f || and imag(u,v) = 0, we chose to rescale the vector x t+1 to 
achieve norm a/II^II ll£ t+1 ll- Thus: 

(5-124) ^ = Mrr) 
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5.5. Robustness to noise and the effect of regularization. We this subsection we make 
a stronger assumption of injectivity: 

Assumption A For every x G H there is a unique X G S 1 ' so that A(X) = A(lx,xJ), 
namely X = [x, xj . 

A heuristic argument (that can be made more precise, but we do not intend to do so here) 
suggests that generically this assumption is satisfied for frames of redundancy 6 or more 
(that is for m > 6n). 

This assumption turns out to be equivalent to a stability bound that we describe next 

Lemma 5.4. The following are equivalent: 

(1) The frame ¥ satisfies hypothesys A. 

(2) ker^n5 2 ' 1 = {0}. 

(3) There is a positive constant A3 > so that 

(5.125) A 3 \\X-Y\\l< ||^(X)-^(F)|| 2 
for allX eS 1 ' 1 andY ES 1 ' . 

Proof 

(l)-v=>(2) is immediate once we notice that S 1 ' — S 1 ' 1 = S 2,1 cf. Lemma 3.6 (5). 

(3)=r*(2) is immediate. 

(2)=>(3). Since S = {W e S 2 ' 1 | {{W^ = 1} is compact, it follows 

(5.126) A3 = M s \\A(W)f = \\A(W )\\ 2 > 

for some W G S 2,1 , W ^ 0. Then since any Z G S 2 ' 1 can be written as Z = WZ^W, with 

W = I ^ZeS, 

\\A(Z)f = \\Z\\l\\A(W)\\ 2 > A 3 \\Z\\l 



Then (|5_125j) follows from noticing that Y-X G S 2 '\ Q.E.D. 

Now we can state the main stability result of the estimators described in this section. 

Theorem 5.5. Fix fi > 2X, x G H zo and let y = a(x) + v . Assume the optimization 
procedure finds a pair u,v G H so that J(u, v, A, /i) < J(x, x, A, /i) . Then 

(5.127) \\{u,vj-{x,x 

Let {u,vj = ctijei, ei] — a 2 [e2,e 2 ], withai,a 2 > and {ei,e 2 } orthonormal set, be its spectral 
decomposition. Assume an oracle provides the global phase e %Vo = rft^Hi ■ Set 




(5.128) x = e vpo y/aiex 

Then 

(5.129) \\x-x\\ 2 < \\{u,v\ - \x,x\\\ x + \a 2 (lu,v] 
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If, additionally, J(u,v, A, /i) < J(0, 0, A,/i) £/jen 

A II II 2 

(5.130) llx-xll 2 < 6— +4||z/|| + ^- 

Proof 

The proof follows from the Lipschitz bound (5.125). Let X = {u,vj and Y = [x, xj. Note 
that 

J 3 (X, A,/i) < J(u,v,X,n) < J(x,x, A,/i) = J 3 ([x,x],A,/i) 
Explicitely this means 

||^4(X - y) - uf + 2A||X|| 1 - 2(2/i - X)a min (X) < \\u\\ 2 + 2A||y|| 1 
Then 
A 3 \\X - Y\\\ < \\A(X - y)|| 2 < (\\A{X - Y) - u\\ + \\v\\f < 2\\A(X -Y)- u\\ 2 + 2||z/|| 2 
Substituting into the previous inequality we obtain 

\\X - Y\\\ ~^\\X- Y\\ x - 4(|M| 2 - (2A - ii)a min (X)) < 

Solving for ||X — Y\\ x we obtain 



2A 

||* - *1i < j- + 2\j ^| + INI" - (2A - /i)a min (X) 

wh ich tog ether with fi > 2A and a miri (X) < for any X G iS 1 ' 1 proves the first inequality 
in (5.127). The second inequality in (refeq: robust) follows from \/a 2 + b 2 < a + b for any 
a,b> 0. 

The second part of the Theorem is obtained as follows. First note [x,x] = ai[ei,ei]. 
Hence 

|p,x] - [x.xllli < ||M - lx,xj\\ x + |a 2 (I«,v])| 

Next we show that ||a; — x\\ < ||[x, x] — [#, x]||i, from where (5.129) follows. Let T = 
{x,xj — {x,xj G S 1 ' 1 . The two (possibly) nonzero eigenvalues are obtained from solving the 
system (3.74). Explicitely we obtain: 

(5.131) a ± = l -tr{T} ± ±y/2tr{T>} - (tr{T})* 

Thus 



\\T\\ t = a + -a_ = v/2tr{T 2 } - {tr{T}) 2 
Direct computations show: 

J. ^r^l^ ll~l|2 II ||2 , frri2^ ll~l|4 , || ||4 ol/~ \|2 

tr{l \— || x 1 1 — \\x\\ , tr{l ) = |p|| + ||x|| — 2\{x,x)\ 

Thus: 

11X11, = (||x|| 2 +||x|| 2 ) 2 -4|(x,x)| 2 
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Recall x has the global phase so that (x, x) > is real and nonnegative. Thus 
((x, x)) 2 and 



KX.X) 






ITHj — 4(x,x) ■ \\x 



x 



which shows 



x 



x\\ < 



\ 1 and thus (5.129 



< \\T\\l 



Finally, (5.130) follows from using the last inequality (5.127) together with a bound for 
1 02] from 

4/x|a 2 | < MX,\,fi) < J 3 (0,A,/i) = \\u\\ 2 . 
Q.E.D. D 

Remark 5. 6. The assumption J(u,v,X, //) < J(0,0, A,/i) is not critical. A different upper 



bound than p. 130) can be obtained from —Afxa2(X) < ^(Y, A,//) 



w\ 



+ 2\\\x\ 



W\ 



\\\x\ 



(5.132) llx-xf < 6— +4||z/|| + ^^ + 

v ' " " _ A 3 " " 4/x 2/i 

which can be made small (to the order of noise level) by choosing a large \x and small X. 

Remark 5.7. Note the results does not require to find the global minimum, but only an 
estimate that brings the criterion J below the level achieved by the true signal. On the other 
hand the result is not unexpected. Similar bounds are known for MLE estimates. 

6. Numerical Analysis 

In this section we present numerical simulations for the Regularized Iterative Least-Square 
algorithm presented in the previous section. 

We generated random frames or redundancy — = 4,6, and 8, as well as complex random 
signals x of size n = 100. All these vectors (frame and signal) are drawn from N(0, 1). Then 
we scale the frame vectors to have norm 1. We set the first component of x to be a positive 
real, and so the global phase becomes uniquely determined. The algorithm stopped when 
A t reached a preset value. In these simulations we choose X mm = 0.01, fi mm = 1 and rate 
7 = 0.8 in n5.11305.114l). 



The magnitude square of signal coefficients a(x) is perturbed additively by Gaussian noise 
with variance a 2 to achieve a fixed signal-to-noise-ratio defined as 



SNR 



ma' 



SNRdB = 10\og w (S NR) [dB] 



We vary SNRdB over 15 values in 5dB increments from -30dB to +40dB. We average 
algorithm performance over 100 noise realizations. 

Figure [l] includes the mean-square error averaged over 100 noise realizations for one fixed 
realization of x, and the C-R lower bound. 

In Figure [2] we plot the bias and variance components of the mean-square error for the 
same results in Figure [Tj Note the bias is relatively small. The bulk of mean-square error is 
due to estimation variance. 

The mean number of iterations varied between 35 and 50, with a higher value for lower 
SNR (-30dB) and a smaller value for higher SNR (40dB). 
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For m = 800 and three values of SNR (-30dB, OdB, and 40dB) Figures (3]§ plot traces of 
a particular noise realization (different for each SNR). In each figure, we plot four character- 
istics: top graph plots the estimation error of the clean signal 101og 10 ||a;' — x|| ; second plot 
presents the smallest eigenvalue of X 1 = [re*, re* -1 ] (which is negative); third plot contains 
the current estimation error of the clean rank-one, that is 101og 10 ||X 4 — [x,x]|| 2 (Frobenius 
norm); bottom graph plots Jo = 101og 10 \\y — ^4(X*)|| . 

We also analyzed the algorithm sensitivity to initialization. If instead of using the eigen- 



pair ( 5.109|5.110 ) we initialize by a random vector x° together with the eigenvalue (5.110). 
Results for the three values of redundancy are presented in Figure [6] for OdB to 40dB range 
oiSNR. 

7. Conclusions 

Novel necessary conditions for complex valued signal reconstruction from magnitudes of 
frame coefficients have been presented. Deterministic stability bounds (Lipschitz constants) 
and stochastic performance bounds (Crameer-Rao lower bound) have been presented. The 
entire analysis has been done canonicaally, that is independent of a pparticular choice of 
basis. Then an optimization algorithm based on the least-square error has been proposed 
and analyzed. The algorithm performance has been compared to the theoretical lower bound 
given by the Cramer-Rao inequality. Remarkably the algorithm performs very well on a 
large range of SNR. In particular, for high SNR, it seems to converge to the correct signal 
every time. This behavior suggests the algorithm presented here is able to track the global 



minimum of (5.105) very well. A future study shall analyze this tracking hypothesis. 
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Figure 1. Mean-Square Error and CRLB bounds for m = 800 (top plot), 
m = 600 (middle plot), and m = 400 (bottom plot) when n = 100. 
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Figure 2 . Bias and variance components of the mean-square error and CRLB 
bounds for m = 800 (top plot), m = 600 (middle plot), and m = 400 (bottom 
plot) when n = 100. 
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Figure 3. Traces for m = 800, n = 100 and SNR = -30dB: 
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Figure 4. Traces for m = 800, n = 100 and SNR = OdB: 
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Figure 5. Traces for m = 800, n = 100 and SNR = 40dB: 
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Figure 6. Mean-Square Error and CRLB bounds for m = 800 (top plot), 
m = 600 (middle plot), and m = 400 (bottom plot) when n = 100 and 
random initialization of x°. 



